library(ggplot2)
  library(reshape2)
  library(dplyr)
  library(tidyr)
  library(GGally)
  library(grid)
  "%&%" = function(a,b) paste(a,b,sep="")
  se <- function(x) sqrt(var(x,na.rm=TRUE)/length(is.na(x)==FALSE))
  source('/Volumes/im-lab/nas40t2/hwheeler/PrediXcan_CV/GTEx_2014-06013_release/transfers/PrediXmod/Paper_plots/multiplot.R')
  my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
  fig.dir <- '~/GitHub/GenArch/GenArchPaper/Figures/'
  rna.dir <- my.dir %&% "gtex-rnaseq/"
  annot.dir <- my.dir %&% "gtex-annot/"
  out.dir <- rna.dir %&% "ind-tissues-RPKM/"
  h2.dir <- my.dir %&% "gtex-h2-estimates/"
tislist <- scan(my.dir %&% '40.tissue.list',sep="\n",what="character")
table1 <- matrix(0,nrow=length(tislist)+1,ncol=6)
ct <- read.table(my.dir %&% 'cross-tissue.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-14.txt',header=T)
n <- ct$N[1]
numexpgenes <- dim(ct)[1]
numexpgenes <- length(ct$local.p[is.na(ct$local.p)==FALSE]) 
meanh2 <- sprintf("%.3f",mean(ct$local.h2,na.rm=TRUE))
semean <- sprintf("%.4f",se(ct$local.h2))
meanandse <- meanh2 %&% " (" %&% semean %&% ")"
pest <-  ct %>%  mutate(pchi=pchisq((local.h2/local.se)^2, df=1, lower.tail=FALSE)) %>% mutate(local.q=p.adjust(pchi,method="BH")) %>% arrange(local.h2) %>% mutate(Qlt05=local.q<0.1)
propsig <- sprintf("%.1f",table(pest$Qlt05)[2]/sum(table(pest$Qlt05,useNA="i"))*100)
numsig <- table(pest$Qlt05)[2]
table1[1,] <- c("Cross-tissue",n,meanh2,propsig,numsig,numexpgenes)
hist(pest$local.p,main="CT")

hist(pest$pchi,main="CT")

for(i in 1:length(tislist)){
  tis <- tislist[i]
  data <- read.table(my.dir %&% 'GTEx.TS.' %&% tis  %&% '.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-14.txt',header=T,sep="\t")  
  n <- data$N[1]
  numexpgenes <- dim(data)[1]
    numexpgenes <- length(data$local.p[is.na(data$local.p)==FALSE]) ##num expressed genes mean(RPKM)>0.1
  meanh2 <- sprintf("%.3f",mean(data$local.h2,na.rm=TRUE))
  semean <- sprintf("%.4f",se(data$local.h2))
  pest <-  data %>%  mutate(pchi=pchisq((local.h2/local.se)^2, df=1, lower.tail=FALSE)) %>% mutate(local.q=p.adjust(pchi,method="BH")) %>% arrange(local.h2) %>% mutate(Qlt05=local.q<0.1) 
  propsig <- sprintf("%.2f",table(pest$Qlt05)[2]/sum(table(pest$Qlt05,useNA="i"))*100)
  numsig <- table(pest$Qlt05)[2]
  meanandse <- meanh2 %&% " (" %&% semean %&% ")"
  tableinfo <- c(tis,n,meanh2,propsig,numsig,numexpgenes)
  table1[i+1,] <- tableinfo
  hist(pest$local.p,main=tis)
  hist(pest$pchi,main=tis)
}

colnames(table1)=c("tissue","n","mean h2","% FDR<0.1","num FDR<0.1","num expressed")
#table1

library(xtable)
tab <- xtable(table1)
print(tab, type="latex",include.rownames=FALSE)
## % latex table generated in R 3.2.0 by xtable 1.8-0 package
## % Wed Aug 10 21:44:04 2016
## \begin{table}[ht]
## \centering
## \begin{tabular}{llllll}
##   \hline
## tissue & n & mean h2 & \% FDR$<$0.1 & num FDR$<$0.1 & num expressed \\ 
##   \hline
## Cross-tissue & 450 & 0.062 & 17.7 & 2995 & 14861 \\ 
##   Adipose - Subcutaneous & 298 & 0.017 & 8.31 & 1408 & 14132 \\ 
##   Adrenal Gland & 126 & 0.037 & 7.84 & 1329 & 14268 \\ 
##   Artery - Aorta & 198 & 0.022 & 8.55 & 1449 & 14634 \\ 
##   Artery - Coronary & 119 & 0.048 & 7.14 & 1211 & 14663 \\ 
##   Artery - Tibial & 285 & 0.022 & 7.44 & 1262 & 14534 \\ 
##   Brain - Anterior cingulate cortex (BA24) & 72 & 0.037 & 9.60 & 1627 & 14904 \\ 
##   Brain - Caudate (basal ganglia) & 100 & 0.042 & 8.19 & 1388 & 14562 \\ 
##   Brain - Cerebellar Hemisphere & 89 & 0.046 & 9.60 & 1627 & 14851 \\ 
##   Brain - Cerebellum & 103 & 0.041 & 9.55 & 1618 & 14675 \\ 
##   Brain - Cortex & 96 & 0.053 & 8.73 & 1480 & 14405 \\ 
##   Brain - Frontal Cortex (BA9) & 92 & 0.037 & 9.26 & 1569 & 14481 \\ 
##   Brain - Hippocampus & 81 & 0.033 & 9.55 & 1618 & 14504 \\ 
##   Brain - Hypothalamus & 81 & 0.020 & 10.85 & 1839 & 14563 \\ 
##   Brain - Nucleus accumbens (basal ganglia) & 93 & 0.040 & 9.35 & 1585 & 14546 \\ 
##   Brain - Putamen (basal ganglia) & 82 & 0.033 & 9.33 & 1581 & 14319 \\ 
##   Breast - Mammary Tissue & 183 & 0.020 & 8.06 & 1367 & 14080 \\ 
##   Cells - EBV-transformed lymphocytes & 115 & 0.045 & 7.80 & 1323 & 14072 \\ 
##   Cells - Transformed fibroblasts & 272 & 0.019 & 7.87 & 1334 & 14193 \\ 
##   Colon - Sigmoid & 124 & 0.024 & 9.31 & 1578 & 14749 \\ 
##   Colon - Transverse & 170 & 0.022 & 9.28 & 1573 & 14217 \\ 
##   Esophagus - Gastroesophageal Junction & 127 & 0.029 & 8.01 & 1358 & 14581 \\ 
##   Esophagus - Mucosa & 242 & 0.026 & 7.20 & 1220 & 14944 \\ 
##   Esophagus - Muscularis & 219 & 0.024 & 7.93 & 1344 & 14805 \\ 
##   Heart - Atrial Appendage & 159 & 0.031 & 8.48 & 1437 & 14401 \\ 
##   Heart - Left Ventricle & 190 & 0.025 & 7.82 & 1325 & 14703 \\ 
##   Liver & 98 & 0.036 & 8.96 & 1519 & 14476 \\ 
##   Lung & 279 & 0.018 & 7.72 & 1308 & 14452 \\ 
##   Muscle - Skeletal & 361 & 0.020 & 7.92 & 1342 & 14382 \\ 
##   Nerve - Tibial & 256 & 0.026 & 7.29 & 1235 & 14347 \\ 
##   Ovary & 85 & 0.043 & 8.98 & 1523 & 13922 \\ 
##   Pancreas & 150 & 0.036 & 8.24 & 1396 & 14678 \\ 
##   Pituitary & 87 & 0.044 & 8.74 & 1481 & 14752 \\ 
##   Skin - Not Sun Exposed (Suprapubic) & 196 & 0.041 & 5.53 & 938 & 15062 \\ 
##   Skin - Sun Exposed (Lower leg) & 303 & 0.027 & 6.02 & 1020 & 14808 \\ 
##   Small Intestine - Terminal Ileum & 77 & 0.046 & 10.22 & 1733 & 13820 \\ 
##   Spleen & 89 & 0.062 & 8.54 & 1448 & 14230 \\ 
##   Stomach & 171 & 0.020 & 9.04 & 1532 & 14259 \\ 
##   Testis & 157 & 0.042 & 9.21 & 1561 & 14779 \\ 
##   Thyroid & 279 & 0.024 & 7.63 & 1293 & 14478 \\ 
##   Whole Blood & 339 & 0.025 & 6.78 & 1150 & 14609 \\ 
##    \hline
## \end{tabular}
## \end{table}